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Abstract 

The Helmholtz equation 

(A + K 2 n 2 )u = f 

with a variable index of refraction n, and a suitable radiation condition at 
infinity serves as a model for a wide variety of wave propagation problems* 
Such problems can be solved numerically by first truncating the given 
unbounded domain and imposing a suitable outgoing radiation condition on an 
artificial boundary and then solving the resulting problem on the bounded 
domain by direct discretization (for example, using a finite element 
method). In practical applications, the mesh size h and the wave number K, 
are not independent but are constrained by the accuracy of the desired 
computation. It will be shown that the number of points per wavelength, 
measured by (Kh)“*, is not sufficient to determine the accuracy of a given 

discretization. For example, the quantity K^h 2 is shown to determine the 

2 

accuracy in the L norm for a second— order discretization method applied to 
several propagation models. 
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Introduction 


The Helmholtz equation 

(1.1) Au + K 2 n 2 u = 0, 

where K is the wave number and n(x) is the index of refraction, describes 
a wide variety of wave propagation phenonraena through an inhomogeneous medium. 
Inhoraogeneities are represented by spatial variations in n(x) and also by 
interfaces and scattering surfaces. Equation (1.1) is fundamental in 
acoustics, in particular, in underwater acoustics ([7], [8]), duct acoustics 
([2], [10], [12]), and acoustical scattering ([5]). In addition, certain 

models of electromagnetic and elastic wave propagation can be described by 
(1.1) ([11], [12]). Vector formulations of (1.1) describe general 

electromagnetic and elastic wave propagation ([15]). Finally, the propagation 
of pulse-like waves can be reduced to an analysis of (1.1) after Fourier 
transforming the time variable ([1]). 

If the wave length X(=2 tt/K) is small relative to the other length 

scales in the problem, solutions to (1.1) can be approximated by asymptotic 
methods. However, if X is of the same order as some characteristic length 
scale, these expansions can break down and the problem must be, in general, 
solved by numerical methods. The methods we are considering are based on 

truncating the domain in which the wave propagation is occurring and imposing 
a suitable outgoing radiation condition on an artificial boundary. The 
resulting problem is then solved on the bounded domain by directly 

discretizing (1.1). Such a method is described in [4], where an efficient 

technique to solve the resulting linear system of equations is also 


introduced 
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In general, radiation conditions do not completely absorb all 
reflections. The total error in' the numerical solution of (1.1) is the sum of 
two errors: the error due to the approximate radiation condition and the 
discretization error due to the approximation of the continuous problem by a 
discrete problem. In this paper we will analyze only the discretization 
errors due to a standard finite element approximation scheme for (1.1), on a 
bounded domain with a suitable radiation condition. 

In any wave propagation problem there are at least three important and 
distinct length scales. These are l - the diameter of the truncated 
computational region, a — the diameter of the region containing the inhomo- 
geneities or other effects which distort free space wave propagation, and the 
mesh size h. Since K has units (length) - ^, this gives three nondimensional 
quantities Ka, K£, and Kh which relate these characteristic lengths to the 
wavelengths. 

(Kh)”^ is the number of grid points per wavelength (up to a constant 
factor) and has been used as a measure of accuracy by many authors (see, for 
example, [2], [7], and [13] and the references contained therein). Ka is 
essentially the number of wavelengths in the inhomogeneous region and is a 
measure of the degree of distortion of the solution from free space wave 
propagation. K£ is a measure of the number of wavelengths in the computa- 
tional domain. It depends on the effectiveness of the radiation boundary 
condition in simulating outgoing radiation and on the positions at which the 
solution is desired. In general, the computational domain is fixed and 
includes all the inhomogeneties. The wave number then varies over some range 
of physical interest. 
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In this paper it will be established that Kh is not a sufficient 
indicator of the truncation error of a discrete approximation to (1.1). The 
arguments, in general, will be given in the context of a finite element 
discretization, nevertheless we expect that similar results are valid for 
finite difference approximations. It will be shown that the discretization 
error depends on both KA and Kh. Thus, when the computational domain is 
fixed, discretization errors will grow as K increases even though the number 
of points per wavelength remains fixed. If a finite element method accurate 
to order m is used, an error bound of o(K m h m will be established for 
errors in the H^-norm. Furthermore, an error bound of 

(1.2) 0(K nH ' 1+C ‘ h m ) 

will be established for errors in the L -norm, where o > 0 depends on both 
the geometry of the problem and the radiation condition. This estimate is 
suboptimal in the sense of approximation theory for the finite element 
subspace. We stress that this analysis is only for the discretization error 
and does not include the errors due to the approximation of the radiation 
condition at a finite boundary. 

Estimate (1.2), with a = 0, was used in [6] in discussing the usefulness 
of the multigrid method to solve the Helmholtz equation. In Section 2, (1.2) 
will be established rigorously in a fairly general setting. It will be shown 
that a = 0 is the most favorable bound and is sharp for a one— dimensional 
model problem but that in general o > 0. The results are obtained from a 
standard finite element error analysis combined with some non-standard lemmas 
bounding the solution in term of the data and K. A reader only interested in 
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the consequences of the theory can skip Section 2 and just read the precise 
statement of Theorem 2.2. Numerical results will be presented in Section 3 
validating the theory in a waveguide geometry* In Section 4 several practical 
consequences of this theory will be discussed* 


2* Error Estimates 

We now outline the theoretical results. We first consider the model 
problem: 

(2.1a) [-A - (K^ + i6K)]u(x) = f(x) x e ft, 

(2.1b) = ° on 3£1, 

where 6 > 0, K > 0, fl is a bounded domain in R^(N = 1, 2, 3) with a smooth 
boundary 3JI and f(x) smooth. 

REMARK 2.1 : The term i5K is introduced so that (2.1) is a well-posed 

boundary value problem. In practical problems this is accomplished by the 

radiation boundary condition. 

To approximate (2.1) we use a finite element method and introduce a 


variational formulation. Let 
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and 


a(u,v) = / [Vu»Vv - (K + i6K)u(x)v(x)]dx 
ft 


(f,v) = / f (x) v(x) dx, 
ft 


then the weak form of (2.1) is 


( 2 . 2 ) 


a(u,v) = (f,v) all v e H X (ft) 


where H (ft) denotes the standard Sobolev space. Given a subspace 

C (ft ) the finite element approximation is the function u* 1 C S* 1 such 
that 


(2.2') a(u h ,v h ) = (f,v h ) for all v h e S h . 

o 

We assume that L functions can be approximated to order h m by elements 
of S h . We can then prove the following theorem. 

THEOREM 2.1: Suppose that u satisfies (2.2) and u e H^^ft). Then 

there exists a unique solution u h of (2.2') provided K 2 h is sufficiently 
small. Furthermore, the following estimates hold for the error e^ 1 = u — u^ 1 

(2.3a) Ie h l < C h m_1 (l + K m )[llul . + y (f)] 

jjl — m v m 

(2.3b) »e h l 2 < C m h m (l + K^^DuH 2 + Y m (f)] 

L L 

where C m depends on m and ft but is independent of K and the data f . 
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Y is given by 
m 


(2.4b) 


Y ( f ) - 
m 


If I 2 
L . 

f 

flf n 

. 2 
H J 

2 + 
K + 1 


+ 1 

■ f « 2 
L + 

t 

" f " j-2 
H J 

K + 1 

j=3 

+ 1 


m even 


m odd 


The sum in (2.4a) ranges over even indices while the sum in (2.4b) ranges over 
odd indices. 

The estimate (2.3b) shows that the L 2 error (normalized by Dull + A m (f)) 

for a scheme of order m grows at least as fast as h m We shall later 

show that in some cases this rate of growth is sharp. For certain classes of 

data f we can also show that y (f) < C fluH 0 where C is independent of 

ra m 

K and f. In these cases we have the estimate 


(2.5) 


ie n l „ < C (l + K™* 1 )h m Hut 


2 m 


and so we have a bound on the relative error lie I' /BuH An example of 

L L 

such a class is data f which can be expanded in a sufficiently rapidly 
convergent series of eigenfunctions of -A in fl. 

A proof of Theorem (2.1) for the model problem (2.1) will be presented 
elsewhere. The proof depends on the finite element analysis of [16] together 
with elliptic estimates and the following bound of the solution in terms of 


the data 
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(2.6) lul 0 < £ Ifl 0 , 

L Z (fl) K L 2 

where C is independent of K and f. For more general problems where 
(2.1b) is replaced by a radiation condition (which can be local or nonlocal, 
see e.g., [3], [8] - [12])), the finite element analysis in [16] has been 

extended to problems with different radiation conditions [9], [10]. However, 

(2.6) is not true, in general, and the strongest bound that we can establish 
is 


(2.7) lul „ < Ifl 

L (£1) K 1_0t L 2 (0) 

where a > 0 depends on the geometry, the dimension of the problem and the 
radiation condition. In such cases we can establish the following bound for 
the error e* 1 


( 2 . 8 ) 


1 e h l < 

I* (ft) 


m 


(K mfl+a + 


l)h m ( lul + Y ( f ) ) • 

i/(ft) m 


A proof of these results will appear elsewhere. We next consider the 
validity of (2.7) for various problems with radiation boundary conditions. It 
can be shown that for the one-dimensional problem 


(2.9) 


“ - K 2 u(x) = f(x) 0 < x < 1 

dx 

u(0) - 0, (l) = iku(l) 


'/ 


where f(x) vanishes near x «* 1, (2.7) holds with a *» 0 
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We next consider the Helmholtz equation in a Cartesian waveguide. Let 
0 = {x e [0,ir], ye and let f = 0 near x = w, and consider the 

problem 


(-A - K 2 )u(x,y) = f(x,y) (x,y) e fl 

( 2 . 10 ) 

u(0,y) = 0, u(x,it) = 0, u (x,0) = 0, u (x.y) = T(u), 

y x 

where T(u) is the global boundary operator for outgoing modes introduced in 
[8] and [10]. In subdomains where f = 0, the solution to (2.9) can be 
expressed as a sum of modes 

» ia .x 

(2.11) u = l q . (y)e 3 

j=0 J 

where 

q^(y) = cos((j+ V 2 )y) 

Oj = / k 2 - (j+ 1/2 ) 2 

For K 2 - (j +V2) 2 > the mode is propagating and outgoing. For 

K 2 - ( 1/2 ) 2 < 0, the j th mode decays exponentially and is called 
evanescent. The values {j+ V 2 } are called cutoff frequencies. When K 
equals a cutoff frequency the solution is not well-posed. 

We can show that (2.7) holds for (2.10) with a * V 2 provided K is 
uniformly bounded away from a cutoff frequency* Furthermore (2*7) holds with 
a = 0 when the solution consists of a finite number of modes. Numerical 
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results for a problem similar to (2.10) will be presented in Section 3. 
Extensions of these results to exterior problems will appear elsewhere. 

2 

For m = 2, the estimate, a = 0, show that as K increases, the L 

o o 

error grows at a rate 0(K h ). In order to show that this growth rate is 
sharp we consider the difference equation in one dimension 

(2.12) u j+l ” ^ u j + u j-l + ^ Uj = 0, 
as a second-order approximation to the equation 

(2.13) u + K 2 u = 0. 

XX 


Equation (2.12) corresponds to discretizing (2.13) with piecewise linear 
elements and lumping the mass matrix (i.e., the terms involving K in the 
bilinear form). It can be seen that the argument below is also valid without 
lumping. 

Solutions to (2.12) are of the form 


where 


. izx. 

=e 1Zjh =e 


= jh, 


(2.14) 


zh = ± Kh[l + 0( (Kh) 2 ) ] . 


If we wish to approximate the outgoing solution (as x +**), the (+) sign 
must be chosen in (2.14) and the approximate solution is 


T 
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la1h iKx (l+0((Kh)^)) 
** e J = e J 


Therefore, the error ej is 


iKx. ix. 0(K 3 h 2 ) 

J [e J - 1], 


e = e - 1 e 


and if we consider a fixed region in x and assume K^h^ small, we obtain 


le.l /lul = 0(K 3 h 2 ). 
3 IT L 


3. Numerical Results 

In order to numerically validate the theory presented in Section 2, we 
consider a model problem 

(3.1) u +u + K 2 u = 0; 0<x<x; 0<y<ir, 

xx yy — — — J — 

with boundary conditions 


u y (x,0) = u(x,ir) = 0 

u x (o»y) = f(y) 

«x<».y> = T < u )» 

where the boundary operator T will be described below. We consider three 
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examples. In example 1, f(y) Is chosen so that the exact solution is 


, . i/ K - .25 x y 

u(x,y) = e cos f , 


and T(u) = i/ K - .25 u. In examples 2 and 3, f is chosen so that the 
exact solution is 


M 


u 


l e i£ j x cos((j+ i)y); t = / K 2 - (j+ 4) 2 , 


j=0 


where M = 4 for example 2 and M = 7 for example 3. The boundary operator 
is the global operator T(u) referred to earlier which was introduced in [8] 

1“ Vi 

for an underwater acoustics propagation model. When is real, the j c 

mode in the solution has no decay in x and is called a propagating mode. 
When Zj is imaginary, the j t ^ 1 mode decays in x and is called evanescent. 

A square N*N grid is used and the equations are solved by the pre- 
conditioned conjugate gradient method described in [4]. Piecewise linear 
elements with lumping are used. Normalized L2 errors for the examples are 
shown in Table I - III for different values of K and N. 

q 2 

In Table I and II the first three entries correspond to K h fixed 
while the first and last two entries correspond to Kh fixed. It is clear 
from the tables that the errors grow almost linearly in K for Kh fixed and 

q 9 

are nearly constant for K J h fixed. In these examples, K is uniformly 
bounded away from the cutoff frequencies and the estimate (2.5) is confirmed 
numerically. (We have observed that this scaling of the error does break down 
as K and N are decreased. This is to be expected from the estimate (2.5) 
as K approaches 0. ) 
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In Table III the first two entries correspond to K 3 h 2 fixed and the 
first, third and fourth entries correspond to Kh fixed. For these entries, 

K is not close to a cutoff frequency and the estimate (2.5) is again 

confirmed. The last three entries contain values of K very near a cutoff 

frequency. In these cases the errors do not scale as predicted and are in 

fact considerably worse. This is because the constant depends on how close 
K is to a cutoff frequency. The errors that would be observed in practice 
depend on the sequence of K values, how close they are to cutoff 

frequencies, and whether the modes close to cutoff are propagating or 
evanescent. 



Table I. 

Results 

for Example 1 


K 

N 

Error 

Kh 

K 3 h 2 

4.16 

65 

.0120 

.204 

.173 

5.45 

97 

.0137 

.178 

.173 

6.60 

129 

.0147 

.162 

.173 

6.24 

97 

.0182 

.204 

.260 

8.32 

129 

.0252 

.204 

.347 



Table II. 

Results 

for Example 2 


K 

N 

Error 

Kh 

K 3 h 2 

4.16 

65 

.0133 

.204 

.173 

5.45 

97 

.0120 

.178 

.173 

6.60 

129 

.0114 

.162 

.173 

6.24 

97 

.0165 

.204 

.260 

8.32 

129 

.0227 

.204 

.347 
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Table III. 

Results 

for Example 3 


K 

N 

Error 

Kh 

K 3 h 2 

4.16 

65 

.013 

.204 

.173 

6.24 

119 

.013 

.166 

.172 

6.24 

97 

.019 

.204 

.260 

8.32 

129 

.025 

.204 

.347 

5.45 

97 

.029 

.178 

.173 

5.55 

97 

.045 

.181 

.183 

6.60 

129 

.036 

.162 

.173 


4. Implications 

We conclude this paper by listing several computational implications of 
the results in Section 2. 

(a) Accuracy evaluations will have to account for the number of 
wavelengths in the computational domain. The number of points per wavelength 
will have to increase with the number of wavelengths to maintain accuracy. 
Thus, the effects of this theory would be expected to become more important as 
new numerical techniques and computer technology make the numerical solution 
to (1.1) feasible for a larger number of wavelengths. For the simple model 
problem (3.1) numerical experiments with a second-order finite difference code 
indicate that we wish to choose the number of points N in each direction to be 
N = .8(Ki) 3/2 to achieve approximately a 7% L 2 accuracy. 

(b) The precise relationship between K and h to maintain a fixed 
accuracy depends on both the order of the discretization scheme, the norm in 
which it is necessary to maintain the accuracy, and also possibly on the 


T 
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geometry and the boundary conditions. The advantages of using higher order 
methods are greater as the number of wavelengths increases. 

(c) Iterative methods for the solution of the linear systems of equations 
obtained by discretizing (1.1) are usually analyzed by studying the 
convergence rate for fixed K as h 0. In practice, K and h are 
constrained by a given accuracy requirement and K increases over some 
interval. Thus, for a second-order method and accuracy determined by the L£ 
norm of the error, these methods should be analyzed for K^h^ fixed (provided 
(2.5) is valid) and K increasing for the propagation models discussed in 
this paper. 
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